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A  CONSTRAINT  ALGORITHM  FOR  MAINTAINING  RIGID  BONDS  IN 
MOLECULAR  DYNAMICS  SIMULATIONS  OF  LARGE  MOLECULES 


I.  INTRODUCTION 

Accurate  integration  of  particle  trajectories  in  molecular  dynamics  (MD)  simu¬ 
lations  requires  timesteps  small  enough  to  resolve  numerical  representations  of  the 
second-order  differential  equation 


(i) 


Common  algorithms  for  integrating  Eq.  (1)  include  Runge-Kutta  and  other  multistep 
methods,  various  predictor-corrector  methods  and  central  differencing  schemes.  The 
Verlet1,  Beeman2,  and  leapfrog  algorithms3’4’5  are  commonly  used  to  integrate  the 
equations  of  motion  in  MD  simulations.  The  simplest  leapfrog  method,  which  inte¬ 
grates  Eq.  (1)  using  only  second-order  central  differences  of  X  and  dX/dt,  has  several 
advantages.  It  is  reversible  with  respect  to  the  independent  variable  (t  in  Eq.  (1)).  It 
requires  fewer  operations  and  less  computer  memory  because  it  does  not  have  to  keep 
track  of  particle  positions  at  several  previous  timesteps.  Finally,  the  leapfrog  algorithm 
tends  to  conserve  energy  better  than  other  methods. 

The  maximum  timestep  size  for  accurately  integrating  the  equations  of  motion  in 
MD  simulations  is  given  by 


St  = 


a 


max 


[»] 


(2) 


where  ma x(dFi/dX)  is  the  maximum  of  the  gradient  of  the  force  of  all  interparticle 
forces  in  the  system,  and  a  is  a  parameter  related  to  the  accuracy  of  the  numerical 
integration.  A  derivation  of  Eq.  (2)  is  given  in  Appendix  1.  In  simulations  of  molecular 

m 

systems,  the  short  timescale  of  intramolecular  forces  usually  determines  the  timestep. 
The  associated  degrees  of  freedom  represent  fast  modes,  such  as  vibrations.  Resolving 
these  modes  is  generadly  not  important  for  simulating  slower  modes  such  as  intermolecu- 
lar  vibrations,  torsional  angle  transitions,  or  molecular  reorientations.  When  the  fastest 
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degrees  of  freedom  can  be  ignored,  the  computational  efficiency  can  be  increased  by 
using  constraint  algorithms  that  eliminate  motions  on  the  fast  scales. 

Various  iterative  and  noniterative  methods  have  been  developed  for  maintaining 
rigid  bonds  in  MD  simulations  of  large  molecules.  For  example,  two  commonly  used 

l 

algorithms  are  SHAKE6”10,  which  is  iterative,  and  the  matrix  method,6  which  is  non¬ 
iterative.  The  matrix  method  inverts  a  matrix  to  solve  for  Lagrangian  multipliers 
that  satisfy  the  constraint  conditions,  and  so  becomes  computationally  expensive  for 
very  large  molecules.  The  SHAKE  algorithm  avoids  matrix  inversion  by  iteratively 
adjusting  particle  coordinates  until  the  system  satisfies  all  the  constraints  to  within  a 
given  tolerance.  In  addition  to  maintaining  rigid  bonds,  constraint  algorithms  must 
counteract  the  increasing  departure  from  the  fixed  distances,  called  constraint  decay, 
resulting  from  the  accumulation  of  numerical  errors.  Iterative  algorithms  counteract 
constraint  decay  implicitly  by  requiring  convergence  to  within  a  specified  tolerance  at 
each  times tep.  Deviations  in  the  constrained  distances  from  their  initial  values  are 
continually  checked  and  corrected.  Noniterative  algorithms  require  an  explicit  scheme 
for  counteracting  constraint  decay  because  there  is  no  inherent  feedback  mechanism  for 
monitoring  changes  in  distance. 

Recently  Edberg  et  al.u  developed  a  noniterative  algorithm  for  maintaining  fixed 
distances  between  particles.  In  conjunction  with  this  algorithm,  they  developed  a 
criterion  to  ensure  that  constrained  distances  only  deviate  from  their  assigned  values 
by  amounts  small  enough  that  the  algorithm  remains  stable.  This  approach  defines 
penalty  functions  that  monitor  constraint  deviations.  When  the  penalty  functions 
reach  some  specified  value,  the  penalties  axe  minimized  by  correcting  the  deviations 
according  to  Gauss’  principle  of  least  constraint12.  By  relaxing  the  constraint  slightly, 
the  computational  cost  is  reduced  because  the  accumulation  of  numerical  errors  is  not 
corrected  every  timestep. 

In  this  paper,  we  present  a  new  algorithm  for  enforcing  holonomic  constraints  that 
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can  be  used  in  conjunction  with  the  leapfrog  algorithm.  We  show  that  it  is  extremely 
efficient,  gives  constraint  forces  explicitly,  and  provides  flexibility  for  computing  con¬ 
straint  forces  according  to  the  most  efficient  programming  structures  for  a  given  prob¬ 
lem  or  computer.  A  Multiple  Constraint  Force  (MCF)  function  is  derived  using  the 
reversible  leapfrog  integration  algorithm  with  constraint  conditions  applied  pairwise 
to  all  restricted  particles.  This  algorithm  is  similar  to  that  developed  by  Memon  et 
al.13  in  that  it  is  based  on  the  leapfrog  integration  algorithm  and  is  iterative.  Rather 
than  solving  for  constaint  forces  in  a  matrix  inversion,  however,  the  MCF  algorithm 
applies  the  two-body  constraint  force  function  iteratively  to  each  restricted  distance 
in  the  polyatomic  molecule.  The  constraint  force  is  added  to  the  other  forces  acting 
on  each  particle,  so  the  value  of  the  constraint  force  function  decreases  with  successive 
iterations.  The  test  simulations  using  this  algorithm  show  that  the  MCF  algorithm 
converges  rapidly,  and  numerical  errors  do  not  accumulate,  so  there  is  no  unstable 
constraint  decay.  Rather,  the  constraint  error  fluctuates  stably.  The  constraint  fluctu¬ 
ations  are  caused  by  small  errors  in  the  constraint  forces  resulting  from:  the  incomplete 
convergence  of  the  constraint-force  functions  when  using  a  small  number  of  iterations; 
approximations  in  the  constraint  force  function  for  the  purpose  of  increasing  computa¬ 
tional  efficiency;  and  numerical  errors  due  to  discretization,  roundoff  and  truncation. 

In  principle,  any  two- particle  constraint  algorithm  could  serve  as  the  basis  of  an 
iterative  multiple-constraint  algorithm  since  the  evaluation  of  any  two-particle  con¬ 
straint  force  can  include  other  constraint  forces  as  external  forces.  For  example,  Singer 
et  al.14  have  developed  an  efficient  and  widely  used  two-particle  constraint  algorithm. 
The  Singer  algorithm,  however,  treats  the  dynamics  of  the  center-of-mass  separately 
from  the  rotational  motion,  thus  treating  the  dynamics  of  the  two  linked  particles  im¬ 
plicitly.  Consequently,  this  method  does  not  permit  a  simple  extension  to  a  system 
with  multiple  constraints. 
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II.  Derivation  of  the  Constraint  Force  Function 

A  constraint  force  function  for  maintaining  rigid  bonds  can  be  derived  by  applying 
the  leapfrog  algorithm  to  solve  Eq.  (1)  numerically.  This  is  done  by  modifying  the 
leapfrog  procedure  to  include  the  forces  of  constraint,  <JF ,  y ,  so  that  the  position  X  and 
the  velocity  V  of  each  particle  are  given  by: 

X^+1  =  X|‘  +  Vtn+i<ft  (3) 


V”+*=V,"-i+(Fr  +  «F”-)it, 

mi 


(4) 


where  £FJy  is  the  constraint  force  for  the  fixed  distance  separating  particles  i  and  j  and 
F”  is  the  sum  of  all  forces  on  particle  i.  The  masses  of  particles  i  and  j  are  given  by 
m,-  and  my,  respectively.  The  superscripts  in  Eqs.  (3)  and  (4)  indicate  that  X  and  V 
are  central  differenced  in  time.  After  n  timesteps  of  size  6t,  V  and  X  sure  computed  at 
times  (n  4-  and  (n  +  l)6t,  respectively. 

The  velocity  of  particle  j  at  step  n  +  £  is 


V"+*  -  V"-* 
j  ~  j 


+  (F?  —  $FJjf) 


6t 


m, 


The  constraint  force  6F?-  must  be  such  that  the  condition 


ix"+i-xr+ii2  =  ^, 


(5) 


(6) 


is  maintained,  where  l0  =  |10|  is  the  fixed  distance  between  particles  i  and  j.  However, 
the  errors  accumulate  at  each  timestep  in  any  real  simulation  so  that 


1x7  -  xn2  =  11, 


(7) 


where  1«  ^  1„.  Equation  6  may  be  written  in  expanded  form  using  the  identities  (3)  to 

(7), 
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(x?  -  x?) + (v;-*  -  ^  ^ )  m2  - 


-*■ 


F? 


MqSFWSt)' 
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where  Mij  =  m,-  +  m;  and  /<  is  incorporated  in  <5F";.  Because  the  constraint  force  acts 
along  the  interparticle  bond,  S has  the  form 

"S  -  ie© 

where  a  is  a  function  of  the  velocities  of  particles  :  and  j,  the  forces  acting  upon 
them,  the  actual  separation  /«  of  the  particles,  and  the  specified  constrained  distance 
!«.  Equation  8  can  be  rewritten  as 


ll=  1«  +  A1  — 


{6t)2a\t 


where 


Letting 


i  i  /  F"  F" 

Ai  =  (v;-‘-v”  *)«+  ~L--1 


mj  m,- 


a(St? 


Eq.  (10)  can  be  rewritten  as  a  quadratic  equation  in  b , 


Solving  for  b ,  an  expression  for  a  follows  from  Eq.  (12): 


a=w[^'-iir±rl+(±:f)L~lA‘)2\  (14) 

Substituting  a  into  Eq.  (9)  (taking  the  negative  root  in  Eq.  (14)),  we  obtain  the  con¬ 
straint  force  function 


«.  -  1  .  k_^l  _  1  ,  (k  • A1)  _  (£2! 

,J  Afy(ft)2  /?  V  ^  1? 


>  < .  .  .  i». 
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A  geometric  interpretation  of  the  constraint  force  function  Eq.  (15)  is  given  in  Appendix 
2.  The  negative  root  in  Eq.  (15)  gives  the  physically  reasonable  value  for  a.  The  square- 
root  term  in  Eq.  (15)  is  the  magnitude  of  the  projection  of  the  displacement  vector  at 
time  t  +  St  onto  the  displacement  vector  at  time  t.  Because  A1  <C  I<,  the  constraint 
force  function  is  given  by  the  approximation 

mimjU  [l  U_A1  III] 

F'3  MijiSty  |_2  +  /2  2/2j’  (16) 

which  for  a  given  timestep  is  less  accurate  them  Eq.  (15),  but  which  is  more  computa¬ 
tionally  efficient  because  it  is  not  necessary  to  compute  a  square  root. 

The  constraint  force  function  <5F y  given  by  Eq.  (15)  is  sufficient  for  maintaining  the 
specified  bond  distance  10  if  5Fy  is  applied  to  a  system  of  two  rigidly  bound  particles. 
For  a  system  of  particles  where  each  particle  may  be  bound  to  more  than  one  other 
particle,  ST y  given  by  Eq.  (15)  is  not  sufficient  for  maintaining  bond  distances  because 
it  assumes  only  a  two-body  constraint  condition  and  neglects  the  cumulative  effect  of 
multiple  connections.  The  constraint  force  function  Eq.  (15)  is,  however,  the  basis  of  an 
iterative  procedure  that  maintains  multiple  rigid  connections  and  also  counteracts  the 
accumulation  of  errors  due  to  discretization  or  approximations  in  the  constraint  force 
functions.  For  a  small  number  of  iterations,  the  computed  constraint  force  function 
may  not  converge  completely  to  its  limiting  value  every  timestep. 

The  resulting  small  error  in  the  computed  value  of  <SFy  combined  with  errors 
resulting  from  discretization  produce  two  types  of  constraint  decay:  deviations  of  the 
bond  lengths  from  their  assigned  values  and  nonzero  velocity  between  constrained  pairs 
of  particles  along  the  direction  of  the  bond.  The  constraint  force  function  Eq.  (15)  has 
two  features  that  counteract  constraint  decay.  These  are  the  conditions  imposed  by 
Eq.  (6)  that  helps  5Fy  to  maintain  an  assigned  bond  length,  and  the  contribution  of 
the  first  term  in  A1  (given  by  Eq.  (11))  to  <5Fy  that  helps  to  minimize  any  relative 
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velocity  components  along  bonds.  These  two  factors  contribute  to  the  stability  of  the 
algorithm  and  give  the  algorithm  the  ability  to  correct  for  the  accumulation  of  small 
errors  resulting  from  discretization  or  any  approximations  in  SFtJ. 

The  second  term  in  Al,  given  by  Eq.  (11),  provides  the  feedback  mechanism  for 
a  convergent  iterative  procedure  for  multiple  connections.  Once  a  constraint  force 
for  a  particular  bond  is  calculated,  it  is  treated  as  an  external  force  acting  on  the 
pair  of  particles.  Thus  at  successive  iterations,  adl  the  constraint  forces  maintaining 
other  bonds  attached  to  either  of  the  two  linked  particles  add  to  the  constraint  force 
maintaining  that  bond. 

The  mathematical  basis  of  the  algorithm  presented  here  is  the  same  as  that  for 
all  other  iterative  constraint  algorithms.  For  each  multiply  connected  particle  in  the 
system,  a  force  AF,-,  the  sum  of  all  constraint  forces  acting  on  the  particle  is  added 
to  the  total  force  F;.  The  quantity  AF,  may  be  given  by 

Ki 

AF  i (17) 

where  the  ary  are  constants  that  specify  the  magnitude  of  each  constraint  force  and 
the  ly  are  vectors  that  specify  the  Kx  bonds  linked  to  each  particle  i.  If  there  are 
N  bonds,  the  constrained  system  has  N  unknowns  ay.  Because  there  are  then  N 
constraint  conditions,  the  system  is  solvable.  Memon  et  al.13  have  discussed  efficient 
matrix  methods  for  determining  these  unknowns  using  constraint  equations  obtained 
using  the  leapfrog  integration  scheme.  Ryckaert  et  al.®  have  pointed  out  that  because 
the  solution  to  the  constraint  system  of  equations  is  unique,  any  convergent  procedure 
that  satisfies  all  geometric  constraint  conditions  by  displacements  of  the  form 


AX.  =  t  ^>0  . 
1=i  mi 


is  equivalent  to  results  obtained  through  direct  solution,  i.e.,  the  matrix  method.  It- 
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erative  constraint  algorithms  differ  in  that  some  adjust  Eq.  (17),  e.g.,  Memon  et  al.13, 
and  others  adjust  Eq.  (IS),  e.g.,  SHAKE6,9.  Each  algorithm,  however,  permits  different 
types  of  programming  structures  which  may  be  optimal  for  particular  types  of  problems 
and  computers.  We  find  that  by  using  an  exact  two-particle  algorithm  as  the  basis  of 
the  MCF  iteration,  the  accuracy  and  convergence  are  both  improved,  and  convenient 
and  efficient  programming  structures  can  be  used  to  evaluate  constraint  forces. 

An  important  aspect  of  the  MCF  algorithm  is  that  it  retains  the  reversibility 
property  of  the  leapfrog  algorithm.  This  is  an  important  guarantee  of  qualitatively  and 
quantitatively  faithful  Hamiltonian  behavior3,4,  as  previously  shown  for  integration  of 
relativistic  charged  particle  orbits  in  electromagnetic  fields. 

Another  important  aspect  of  the  MCF  algorithm  is  that 


(5F;? 


for  increasing  (i),  independent  of  the  order  in  which  the  <5FtJ  are  computed.  The  su¬ 
perscript  (i)  in  Eq.  (19)  designates  the  value  of  the  constraint  force  function  at  the 
ith  iteration.  This  property  provides  flexibility  for  computing  6Ft]  at  each  iteration 
because  it  is  possible  to  partition  the  calculation  of  constraint  force  functions  accord¬ 
ing  to  what  is  more  efficient  for  a  given  problem  or  computer.  Further,  F|;j)  can  be 
evaluated  unambiguously  before  the  iteration  begins. 

The  constraint  force  function  F,;  is  an  explicit  function  of  all  forces  acting  on  the 
linked  particles.  An  important  optimization  permitted  by  the  MCF  algorithm,  as  well 
as  other  iterative  constraint  procedures  such  as  SHAKE6,9,  is  to  set  the  forces  between 
all  linked  particles  to  zero.  This  step  reduces  the  magnitude  of  the  maximum  gradient 
of  the  constraint  forces  evaluated  for  the  system  and  larger  timesteps  can  be  used. 
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III.  Procedure  for  Calculating  Constraint  Forces 

A  flow  chart  describing  the  MCF  algorithm  is  displayed  in  Fig.  1.  For  each  particle 

in  the  system: 

(1)  Compute  the  sum  of  all  the  forces  acting  on  each  particle  due  to  the  other  particles 
in  the  system. 

One  important  feature  of  this  stage  of  the  calculation  is  that  the  force  between 
constrained  particles  is  set  to  zero.  This  is  not  necessary  in  principle  because  the 
constraint  force  counteracts  any  mutual  interaction.  However,  short-range  interactions, 
typically  at  bond-length  separations,  are  relatively  large  and  result  in  a  large  value  for 
terms  containing  1<  •  Al«  in  the  constraint  force  function.  Because  the  gradient  of  the 
force  is  large  at  short  range,  if  we  included  these  short-range  interactions,  we  would 
need  a  very  small  timestep. 

(2)  Compute  the  displacement  vector  A1  given  by  Eq.  (11)  for  each  bond  in  the  system. 

The  computation  of  A1  for  the  various  bonds  must  be  independent  of  the  order 
of  computation  in  any  convergent  iterative  scheme.  That  is,  the  convergence  does 
not  depend  on  whether  A1  is  computed  for  the  different  bonds  according  to  a  specific 
procedure.  Although  different  modes  of  computing  the  displacement  vector  A1  might 
result  in  different  sets  of  values  of  A1  for  the  first  iteration,  successive  iterations  should 
diminish  this  difference. 

(3)  The  constraint  force  for  maintaining  each  constrained  distance  is  computed  ac¬ 
cording  to  the  constraint  force  function  Eq.  (15)  or  Eq.  (16). 

(4)  The  constraint  force  for  each  constrained  distance  is  added,  with  the  appropriate 
sign,  to  the  total  force  acting  on  the  two  constrained  particles. 


(5)  If  the  maximum  constraint  force  calculated  in  the  ith  iteration  is  less  than  the  con¬ 
vergence  condition,  or  if  the  iteration  counter  exceeds  the  iteration  limit,  proceed 
to  integrate  the  equations  of  motion;  otherwise,  repeat  steps  (2)  to  (4). 


The  development  of  the  constraint  forces  during  the  first  two  iterations- is  illus¬ 
trated  in  Fig.  2.  For  this  example,  we  have  chosen  to  compute  the  displacement  vector 
sequentially  from  top  to  bottom.  As  the  timestep  St  decreases,  the  order  of  evaluation 
has  less  of  an  effect  on  the  total  force  eventually  applied  to  each  particle.  The  constraint 
conditions  are  coupled  to  one  another  by  treating  constraint  forces  as  external  forces 
for  subsequent  calculations.  Thus,  while  <SF^  depends  only  on  external  forces, 


depends  on  external  forces  and  the  constraint  force  ST^' .  On  the  second  iteration,  F 
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depends  on  the  external  forces  and  F^.  This  coupling  leads  to  overall  convergence  of 
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IV.  Tests  of  the  Algorithm 

We  conducted  a  series  of  simulations  to  test  our  algorithm.  The  purpose  of  these 

simulations  was  to 

(1)  examine  the  stability  and  accuracy  of  the  algorithm; 

(2)  test  those  aspects  of  the  algorithm  that  contribute  to  its  general  efficiency; 

(3)  generate  an  initial  version  of  a  computer  program  based  on  this  algorithm  to 

be  used  for  large  scale  simulations. 

The  model  consisted  of  a  12  x  12  array  of  rigid  tetra-atomic  molecules  in  three- 
dimensional  space.  The  configuration  of  each  tetra-atomic  molecule  is  shown  in  Fig.  3. 
The  bond  lengths  and  angles  were  fixed  at  the  normal  carbon-carbon  single  bond  length 
of  1.54  A  and  at  the  tetrahedral  angle,  respectively.  Non-bonded  interactions  were 
given  by  a  Leonard- Jones  potential  with  parameters  taken  from  atomic  nitrogen16, 
e  =  0.5143  x  10“14  erg  and  <t  =  3.310  A.  The  boundary  conditions  were  periodic  in 
the  plane  of  the  initial  array.  The  system  was  confined  in  the  third  dimension  by  a 
pair  of  walls  with  the  same  LJ  parameters  as  the  particles  themselves.  Each  molecule 
in  the  system  was  given  a  random  initial  velocity,  and  the  motion  was  calculated  with 
times teps  varying  between  2  x  10~ia  and  2  x  10“us.  The  deviation  of  actual  distances 
(/«)  from  their  proper  values  (/„)  was  monitored  during  each  calculation.  The  number  of 
distances  deviating  by  more  than  1%  and  the  identity  of  the  deviating  pairs  of  particles 
were  stored. 

Figure  4  shows  the  system  in  its  initial  configuration,  i.e.,  t  —  0.  Figuxes  4a,  4b 
and  4c  are  the  XZ,  YZ,  and  XY  projections  of  the  system,  respectively,  at  t  =  0.  The 
“periodic  boundaries"  are  at  0  A  and  60  A  along  X  and  Y,  while  the  wadis  in  Z  are  at 
-4  A  and  +8  A. 

The  energy  and  constraint  deviations  are  displayed  in  Fig.  5a  and  Fig.  5b,  respec¬ 
tively,  for  a  calculation  of  42000  steps  of  5  x  10-I6s  each,  using  the  complete  constraint 
force  function  Eq.  (15).  The  constraint  force  was  evaluated  with  five  iterations  at  each 


timestep.  The  mean  velocity  of  the  atoms  and  centers  of  mass  of  the  molecules,  along 
with  the  corresponding  standard  deviations,  are  given  in  Table  1  (state  1).  Also  shown 
are  the  values  of  the  mean  and  standard  deviation  of  the  total,  kinetic,  and  potential 
energies.  The  reader  will  note  the  stability  and  rather  narrow  range  of  fluctuations  in 
all  these  values.  Figure  6b  shows  the  number  of  constrained  bonds  that  deviate  from 
their  proper  values  by  more  than  1%.  The  maximum  deviation  found  under  these  con¬ 
ditions  was  about  2%.  The  average  number  of  deviations  per  timestep  is  about  3,  and 
the  maximum  recorded  is  12,  out  of  720  bonds  in  the  whole  system.  The  identity  of  the 
constrained  particles  responsible  for  these  deviations  does  not  stay  the  same  for  many 
timesteps.  Both  the  small  number  and  the  fluctuating  identity  of  deviations  demon¬ 
strate  the  ability  of  the  MCF  algorithm  to  counteract  constraint  decay.  Under  these 
conditions,  we  prefer  to  speak  of  “constraint  fluctuation,”  since  there  is  no  tendency 
for  the  number  or  magnitude  of  the  deviations  to  grow. 

A  similar  calculation  of  20000  steps  using  the  approximate  version  of  the  constraint 
force  function  given  by  Eq.  (16)  is  presented  in  Fig.  6.  The  state  of  the  system  (state 
1,  in  Table  1)  remained  the  same,  well  within  the  standard  deviations.  The  number  of 
constraint  violations  at  the  1%  level  is  a  factor  of  thirty  larger,  but  the  magnitude  of 
any  single  deviation  is  about  the  same,  the  largest  deviation  observed  being  about  2%. 
In  this  case  too,  the  identity  of  the  deviating  bonds  persisted  over  approximately  100 
timesteps. 

Figure  7  shows  the  energies  of  a  system  described  by  the  parameters  listed  for  state 
2  in  Table  1  over  3000  steps  of  1  x  10“15s  (3  x  10-12s  total  time)  with  the  constraint 
force  function  evaluated  over  10  iterations  in  reverse  order  from  the  other  runs  shown 
thus  far.  This  shows  that  the  stability  of  the  algorithm  does  not  depend  on  the  order  of 
the  evaluation  of  the  constraint  force.  A  comparison  of  the  calculation  shown  in  Fig.  7 
and  a  calculation  starting  from  the  same  initial  state,  but  with  the  constraint  forces 
calculated  in  the  standard  sequence,  is  shown  in  Table  2.  Although  the  centers  of  mass 
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of  the  molecules  move  4.30  A  on  average,  the  final  states  corresponding  to  the  “forward” 
and  “backward”evaluation  differ  by  1  part  in  105  in  the  mean  displacements.  The  mean 
velocities  and  energies  of  the  particles  and  molecules  show  a  similar  agreement.  Table 
2  also  presents  the  results  for  forward  versus  backward  runs  at  timesteps  of  1  x  10~14s 
and  .2-5  x  10~16s.  Although  the  total  energy  is  conserved  for  all  runs,  the  1  x  10-14s 
times tep  is  too  large  under  these  circumstances  and  there  is  a  significant  dependence 
on  the  order  of  evaluating  the  constraint  force  functions. 

Figure  8  shows  the  absolute  magnitude  of  the  constraint  force  as  a  function  of 
iteration  for  one  timestep  of  2  x  10“14  sec,  in  state  2  of  Table  1.  The  solid  lines  show 
the  values  of  the  largest  constraint  forces  in  the  entire  system,  and  the  dashed  line 
shows  the  average  constraint  force.  The  average  constraint  force  is  approximately  an 
order  of  magnitude  less  than  the  maximum  constraint  force.  Also,  both  the  average  and 
the  maximum  fall  off  exponentially  at  about  the  same  rate  of  two  orders  of  magnitude 
in  ten  iterations.  Furthermore,  the  standard  deviation  about  the  mean  constraint  force 
is  of  the  same  order  as  the  mean  force  itself  (Table  3).  Thus  max  ( 6F, appears  to 
be  a  good  diagnostic  for  the  convergence  of  the  constraint  force  evaluation.  Figure 
8  also  shows  that  the  constraint  force  does  not  decrease  monotonically  with  iteration 
number,  and  that  the  identity  of  the  pair  of  particles  requiring  the  maximum  constraint 
force  changes  with  the  number  of  iterations.  An  examination  of  the  behavior  of  the 
constraint  force  function  over  many  timesteps  shows  that  those  particle  pairs  requiring 
the  largest  constraint  forces  at  early  iterations  change  identity  slowly,  over  several 
hundred  timesteps  under  the  present  conditions  (state  2),  and  that  the  magnitude 
of  these  forces  also  changes  slowly.  These  last  observations  suggest  that  a  predictor- 
corrector  version  of  the  algorithm,  in  which  a  previous  evaluation  of  the  constraint  force 
is  used  as  the  initial  estimate,  would  be  very  efficient.  Preliminary  tests  support  this 
idea. 
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V.  Conclusion 

We  have  derived  and  tested  a  new  algorithm  for  constrained  dynamics.  The  nu¬ 
merical  basis  of  the  algorithm  is  established  by  applying  a  constraint  condition  to  the 
equations  of  motion  in  the  central-difference  leapfrog  integration  scheme,  and  solving 
for  the  forces  of  constraint.  A  geometric  interpretation  is  provided  in  Appendix  2  that 
provides  physical  insight  to  the  terms  of  the  constraint  force  function.  An  approximate 
version  of  the  constraint  force  function  that  eliminates  the  need  to  evaluate  a  square 
root  saves  approximately  10%  in  execution  time.  However,  the  ultimate  cost  of  the 
decrease  in  accuracy  has  yet  to  be  assessed.  The  algorithm,  which  is  exact  for  two- 
body  systems  such  as  rigid  rotors,  is  applied  iteratively  in  simulations  of  polyatomic 
molecules.  At  the  end  of  each  iteration,  the  calculated  constraint  forces  are  added  to 
the  forces  acting  on  the  system.  The  largest  constraint  force  falls  off  exponentially  with 
iteration  number,  and  may  be  used  as  a  diagnostic  of  convergence. 

We  are  developing  and  testing  several  simple  extensions  and  improvements  of  the 
MCF  algorithm  that. seem  worthy  of  consideration  for  future  implementation.  Storing 
the  composite  constants  {ary}  for  each  of  the  constrained  distances,  rather  than  the 
total  constraint  force  derived  from  these  constants,  allows  calculation  of  better  approx¬ 
imations  to  the  constraint  force  in  the  first  iteration  at  the  next  timestep,  when  the 
geometry  has  changed  somewhat.  Using  the  total  constraint  force  suffers  from  the  fact 
that  the  orientation  of  the  atoms  can  be  different  and  hence  the  vector  constraint  force 
must  change,  even  though  its  magnitude  along  the  line  of  centers  may  be  identical. 
Using  the  constants  {a,-/}  has  the  added  advantage  that  only  one  scalar  constant  has 
to  be  stored  for  each  constrained  distance  rather  than  the  three  components  of  the 
corresponding  vector  force. 

If  the  constraint  constants  {ay}  are  kept  for  the  two  previous  steps  rather  than  one, 
a  simple  extrapolation  to  the  next  timestep  should  give  a  good  (linear)  approximation 
to  the  inevitable  changes  in  the  values  of  {ay}  that  occur.  This  extrapolation  should 
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work  best  for  the  slowly  changing  parts  of  {&,_,}  for  which  the  most  iterations  are 
required.  The  iteration  procedure  treats  the  rapid  changes  in  the  constraint  constants 
that  arise  from  close  interactions  and  impulses  with  inherent  efficiency.  Thus  it  should 
be  possible  to  reduce  the  number  of  iterations  required  appreciably,  with  improved 
accuracy,  at  essentially  no  additional  computer  cost. 

The  approximation  to  the  square  root  used  in  Eq.  (16)  reduces  the  overall  com¬ 
putation  noticeably;  but,  as  the  differences  between  figures  5  and  6  imply,  it  is  rather 
crude.  Since  a  few  iterations  are  necessary  to  obtain  the  constraint  force  anyway,  the 
solution  can  be  obtained  to  the  roundoff  limit  by  using  a  quadratically  convergent  iter¬ 
ative  approximation  to  the  square  root.  This  procedure  should  be  more  efficient  than 
Eq.  (15)  and  equally  accurate. 

The  constrained  equations  of  motion  obtained  here  are  ultimately  the  same  as 
those  obtained  by  other  constraint  algorithms,  an  especially  clear  summary  of  which 
is  presented  by  Levitt  and  Meirovitch17.  A  constraint  force  evaluated  at  a  previous 
timestep  can  be  stored  as  a  particle  attribute  to  be  used  for  subsequent  information 
processing,  as  in  predictor-corrector  schemes18  or  interparticle  state  transitions19. 

Although  much  of  our  approach  has  been  influenced  by  the  work  of  Edberg  et 
al.n,  particularly  the  concept  of  penalty  functions,  such  functions  are  not  necessary  in 
our  formulation.  We  attribute  this  to  the  remarkable  stability,  demonstrated  in  part 
IV  of  this  paper,  of  the  reversible  leapfrog  integration  scheme.  We  note,  however,  that 
energy  conservation  is  not  a  sufficient  indicator  of  the  accuracy  of  the  algorithm.  A 
better  test  is  that  the  statistical  properties  sought  should  not  be  unduly  perturbed 
by  the  choice  of  timestep  size.  We  have  also  derived  a  simple  relation  between  the 
maximum  force  gradient  experienced  by  the  system,  and  the  maximum  timestep  to 
integrate  the  equations  of  motion  accurately  (see  Appendix  1). 

Finally,  we  have  written  a  computer  program  that  uses  the  MCF  algorithm  to 
conduct  molecular  dynamics  simulations  of  large  assemblies  of  molecules.  Among  the 
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issues  we  are  addressing  are:  development  of  optimum  programming  structures  based 
on  the  two-paxticle  constraint  force  approach  of  the  MCF  algorithm;  development  of 
efficient  vector  procedures  based  on  the  property  that  the  MCF  method  converges 
independent  of  the  order  in  which  the  Fij  axe  computed;  and  extending  the  MCF 
method  to  predictor-corrector  schemes  for  computing  constraint  forces  as  described 
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Figure  1.  Sequence  of  operations  required  for  maintaining  constrained  distances 

Eqs.  (11)  and  (15). 


FIRST  ITERATION  SECOND  ITERATION 


constraint  forces.  The  quantity  <5F,y  is  the  constraint  force  for  maintaining  bond  ( ij ) 
at  the  ith  iteration. 
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TIMESTEPS  (5  x  10-16s/step) 

Figure  5a.  Energy  as  a  function  of  time.  Upper  curve  -  total  energy  (kinetic 
tential).  Lower  curve  -  total  kinetic  energy.  (The  constraint  force  is  given  bv 
expression  Eq.  (15)). 
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Figure  6a.  Energy  as  a  function  of  time.  Upper  curve  -  total  energy  (kinetic  +  poten¬ 
tial).  Lower  curve  -  total  kinetic  energy.  (The  constraint  force  is  given  by  approxoimate 
expression  Eq.  (16)). 
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TIMESTEPS  (5x  10"16s/step) 

Figure  6b.  Number  of  constraint  deviations  as  a  function  of  time.  (The  constraint  force 
is  given  by  approximate  expression  Eq.  (16)). 
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Figure  7.  Total,  kinetic,  and  potential  energy  of  system  as  a  function  of  time  corre¬ 
sponding  to  a  reversal  of  the  order  of  evaluating  the  constraint  forces.  (Exact  expression 
Eq.  (15),  10  iterations). 
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Appendix  1. 


Derivation  of  Eq.  (2):  Maximum  Timestep  for  Accurately 
Integrating  Equations  of  Motion 


The  trajectory  of  a  particle  of  mass  m  in  a  force  field  F(X)  is  determined  by 


(A1  -  1) 


where  X  is  the  position.  In  general,  F(X)  is  a  nonlinear  function  of  X.  However,  an 
accurate  approximation  of  Eq.  (Al-1)  by  difference  equations  depends  on  the  form  of 
F(X)  in  some  small  neighborhood  SX.  Applying  a  perturbation  to  Eq.  (A  1*1)  and 
expanding  F  in  a  Taylor  series  about  X ,  we  obtain 


rn~X£FL-)'  =  F(*  +  SX)  »  F(x)  +  Sx%- 


(.41  -  2) 


Subtracting  Eq.  (Al-1)  from  Eq.  (Al>2),  we  obtain 


(.41  -  3) 


The  quantity  d?{8X)Jdt2  can  be  approximated  by  applying  central  differences  to  <5 X, 


cP(SX)  ^  ^ SX(t  +  St)  -  2 SX(t)  4-  SX(t  -  St) 


(-41  -  4) 


for  St  sufficiently  small.  Bounds  on  the  size  of  St  are  determined  by  expanding  6X(t±St) 


in  a  Taylor  series  about  t, 


dX  ( St )2  d2(SX) 

sx(t  ±  st) » sx(t)  ±  st-£-  +  1  ;  (  J 


dt  2!  dt2 


(.41-5) 


{sty  d3{sx)  (sty  d\6x) 

3!  dt3  4!  dt 4 


Substituting  Eq.  (Al-5)  into  the  right  side  of  Eq.  (Al-4),  we  obtain 


SX(t  +  St)  -  2 SX(t)  +  SX(t  -St)  _  d?(6X)  St 2  d4(SX) 

St 2  %  dt 2  +  12  dt4 


(.41-6) 


Equation  Al-6  shows  that  the  approximation  given  by  Eq.  (Al-4)  is  accurate  only  if 


St2  d4(SX)  /d*(*A')\ 

12  dt 4  V  ^  )  ~7, 


(.41  -  7) 


where  7  C  1.  It  follows  from  Eq.  (Al-3)  that 


cfijSX)  _  1  <P(6X)  dF 
dt 4  ~  m  dt 2  dX' 


(.41-8) 

Substituting  Eq.  (Al-S)  into  Eq.  (Al-7)  and  letting  a  =  (12m7)1,/2,  we  obtain 


St  = 


i£ 


r 


(.41-9) 


which  is  Eq.  (2)  in  the  text. 
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Appendix  2. 


Geometric  Interpretation  of  Eq.  (15) 


The  constraint  force  maintaining  a  fixed  bond  length  between  two  particles  given 
by  Eq.  (15)  can  also  be  obtained  by  considering  the  relative  motion  of  these  particles. 
The  motion  of  a  particle  of  mass  m*  relative  to  a  particle  of  mass  m;  is  equivalent 
to  the  motion  of  a  particle  of  mass  n  =  ,  -f  m;)  moving  relative  to  a  fixed 

point  in  space  in  the  center-of-mass  coordinate  system.  We  therefore  consider  the  force 
required  to  constrain  the  motion  of  a  particle  of  mass  n  to  a  fixed  orbit  about  a  point 
in  space. 

The  constraint  force  function  Eq.  (15)  may  be  expressed  as  the  sum  of  two  forces. 


6F  ij  —  6  Fci  +  SF  c2  ■ 


(A2-  1) 


The  force  £Fei  counteracts  the  external  forces  that  tend  to  cause  deviations  from  the 
fixed  bond  length.  The  force  <5FC2  both  gives  the  bound  particle  a  trajectory  consistent 
with  the  rigid-bond  constraint  and  counteracts  the  influence  of  accumulated  numerical 
errors  that  cause  the  bond  length  to  deviate  from  its  fixed  value.  The  force  <5FC] ,  given 
in  terms  of  the  geometric  quantities  defined  in  Section  III, 


*Fcl 


A*  U  •  A1  / l«  \ 

( sty  i<  \i')' 


(.42-2) 


is  explained  by  Fig.  Al-1.  The  quantity  Al,  given  by  Eq.  (11),  is  the  displacement 
vector  of  the  particle  in  the  center-of-mass  system.  As  can  be  seen  in  Fig.  A2-1,  the 
quantity  1  sc  is  the  position  of  the  particle  at  time  t +  <5t  if  the  constraint  force  Eq.  (.42- 
2)  is  not  applied.  The  force  6FC2,  where 


6  Fc2  = 


(sty 


-0<-ic), 


(-42  -  3) 
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is  also  explained  by  Fig.  A 1-1.  Note  that  in  principle  |1{  |  =  |10|,  but  in  actual  simulations 
the  accumulation  of  numerical  errors  destroys  this  equality.  The  quantity  lc  is  the 
projection  of  the  displacement  vector  at  time  t+6t  onto  the  displacement  vector  at  time 
t.  Because  we  desire  a  fixed-orbit  trajectory  of  orbital  radius  |I0|,  the  force  Eq.  (A2-3) 
must  be  applied  with 


u  =  U 


ll  +  ~  (A')2- 


(.42-4) 


Figure  Al-1  gives  a  geometric  description  of  Eq.  (A2-4).  Note  that  in  addition  to 
adjusting  the  trajectory  in  accordance  with  the  constraint  distance,  the  constraint 
force  also  adjusts  the  trajectory  to  correct  for  any  deviations  due  to  numerical  error. 
Substituting  Eq.  (A2-4)  into  Eq.  (A2-3)  and  adding  Eq.  (A2-2),  we  obtain  Eq.  (15). 
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Table  1.  Average  states  of  system  used  to  test  stability  and  accuracy  of  algorithm. 
State  1  refers  to  figures  6  and  7.  State  2  refers  to  figures  8  and  9.  Velocities  are 
averaged  over  the  entire  system  at  one  times tep.  Energies  are  averaged  over  an  entire 
run.  State  1  energies  averages  are  over  20000  timesteps;  state  2  energy  averages  are 
over  3000  timesteps. 


Mean  Particle 
Velocity 
(cm/s) 

Mean  Molecular 
Velocity 
(cm/s) 

Total 

Energy 

(erg? 

Kinetic 

Energy 

(erg) 

Potential 

Energy 

(erg) 

State  1:  1.98x10? 

±8.88xio4 

1.27x10? 

±4.91xl04 

3.66xl0-10 

3.32x10*1° 

3.38x10*1° 

State  2:  2.05X104 

±9.41x103 

137x10* 

+5.52X103 

-3.65x10*12 

3.4 1x10* 12 

-7.05x10*  I2 

Table  2.  Average  displacements  of  particles,  ( 6rat ),  and  centers-of-mass,  ( 6rcm ).  be¬ 
tween  forward  (F)  and  reverse  (R)  evaluation  of  the  constraint  force  functions.  The 
total  period  simulated  was  3  x  10” 12  sec  in  ail  cases. 


timestep' 

1x10* 

14  s 

1x10- 

■15  s 

1x10' 

direction 

F 

R 

F 

R 

F 

<irat> 

7.60±14.92 

7.45±14.96 

5.38±12.96 

5.38+12.96 

5.98+14.07 

5.77±  7.49 

5.67+  7.47 

4J1±,  6.52 

431+  6.52 

4.44+  6.71 

16  s 

R 

5.98+14.07 
4.44+  6.71 
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Table  3.  Maximum  and  average  constraint  force 


a  function  of  iteration. 


ration 

Bond  Index 

Maximum  <5Fi; 

Average  <5Fi; 

1 

269 

7.049  x  10-6 

9.735  x  10“' 

2 

212 

2.674  x  10_<5 

4.3S7  x  10~7 

3 

267 

1.622  x  10~6 

2.690  x  10“7 

4 

149 

1.044  x  10~6 

1.605  x  10~7 

5 

149 

6. 896  x  10"7 

9.524  x  10“8 

6 

149 

4.4S9  x  10-7 

5.692  x  10“8 

t 

149 

2.S79  x  10“7 

3.397  x  10~8 

8 

149 

1.S25  x  lO"7 

2.02S  x  10~8 

9 

149 

1.145  x  10“7 

1.205  x  10~8 

10 

149 

7.126  x  10“8 

7.120  x  lO"9 

11 

149 

4.400  x  10"8 

4.201  x  10~9 

12 

149 

2.699  x  10“8 

2.473  x  lO"9 

13 

149 

1.645  x  10~8 

1.454  x  10~9 

14 

241 

1.031  x  lO"8 

8.555  x  10~10 

15 

241 

6.799  x  10~9 

5.030  x  10~10 

16 

241 

4.4S6  x  10-9 

2.964  x  10_1° 

17 

241 

2.960  x  10-9 

1.753  x  10~10 

IS 

241 

1.953  x  lO"9 

1.040  x  10_1° 

19 

241 

1.2S9  x  lO"9 

6.175  x  10~n 

20 

241 

S.503  x  10-1° 

3.67S  x  10"11 
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